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Abstract 



Wc present a method to solve the Hehiiholtz equation for a non-homogeneous mem- 
brane with Dirichlet boundary conditions at the border of arbitrary two-dimensional 
domains. The method uses a collocation approach based on a set of localized func- 
tions, called "little sine functions", which are used to discrctize two-dimensional 
regions. Wc have performed extensive numerical tests and wc have compared the 
results obtained with the present method with the ones available from the literature. 
Our results show that the present method is very accurate and that its implemen- 
tation for general problems is straightforward. 
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1 Introduction 



This paper focuses on the solution of the inhomogeneous Helmholtz equation 



over an arbitrary two dimensional membrane B, with Dirichlet boundary con- 
ditions at the border, dB. u{x,y) is the transverse displacement and E — cu^, 
u being the frequency of vibration of the membrane. 
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This problem has been considered in the past by several authors, using dif- 
ferent techniques: for example, Masad in [T] has studied the vibrations of a 
rectangular membrane with linearly varying density using a finite difference 
scheme and an approach based on perturbation theory; the same problem 
was also considered by Laura and collaborators in |2j, using an optimized 
Galerkin-Kantorovich approach and by Ho and Chen, [3], who have used a 
hybrid method. Recently Reutskiy has put forward in |1] a new numerical 
technique to study the vibrations of inhomogeneous membranes, the method 
of external and internal excitation. Finally, Filipich and Rosales have studied 
the vibrations of membranes with a discontinous density profile. 

In this paper we describe a different approach to this problem and compare its 
performance with that of the methods mentioned above. Our method is based 
on a collocation approach (see for example [7]) which uses a particular set 
of functions, the Little Sine funetions (LSF), introduced in |8|9] . to obtain a 
discretization of a finite region of the two-dimensional plane. These functions 
have been used with success in the numerical solution of the Schrodinger 
equation in one dimension, both for problems restricted to finite intervals 
and for problems on the real line. In particular it has been observed that 
exponential convergence to the exact solution can be reached when variational 
considerations are made (see [8|9] ). In ref. [9], the LSF were used to obtain a 
new representation for non-local operators on a grid and thus numerically 
solve the relativistic Schrodinger equation. An alternative representation for 
the quantum mechanical path integral was also given in terms of the LSF. 

Although Ref.jH] contains a detailed discussion of the LSF, I will briefly review 
here the main properties, which will be useful in the paper. Throughout the 
paper I will follow the notation of [8]. 

A Little Sine Function is obtained as an approximate representation of the 
Dirac delta function in terms of the wave functions of a particle in a box 
(being 2L the size of the box). Straightforward algebra leads to the expression 



s^fh N x)^— I sin((2Ar + l) x-{x)) _ cos ((2Ar + l)x+(x)) 
' ' 2A^ I sinx_(x) cosx+(a^) 

where x±{x) = 2Wh^^ -'- ^^)- '^^^ index k takes the integer values between 
~N/2 + 1 and A^/2 — 1 (A^ being an even integer). The LSF corresponding to 
a specific value of k is peaked at Xk = 2Lk/N = kh, h being the grid spacing 
and 2L the total extension of the interval where the function is defined. By 
direct inspection of eq. (|2| it is found that Sk{h, N, Xj) = 5kj, showing that 
the LSF takes its maximum value at the fc*^ grid point and vanishes on the 
remaining points of the grid. 

It can be easily proved that the different LSF corresponding to the same set 
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are orthogonal [8]: 



L 

j Sk{h, N, x)sj{h, N, x)dx = h 5kj (3) 

-L 

and that a function defined on a; G {—L, L) may be approximated as 



N/2-1 

f{x)^ fixk) Skih,N,x) . (4) 

k=-N/2+l 



This formula can be applied to obtain a representation of the derivative of a 
LSF in terms of the set of LSF as: 



dsk{h,N,x) dskih,N,x) 



dx dx 

1 



Sj {h,N,x) = J2 4i Sj{h, N, X ) 



d^Skih.N.x) y-^ d^Skih.N.x) 



dx'^ ~^ dx"^ 



kj V""' ! 5 

j 



(5) 



where the expressions for the coefficients c^^j can be found in [S]. Although 
eq. Q is approximate and the LSF strictly speaking do not form a basis, 
the error made with this approximation decreases with N and tends to zero 
as tends to infinity, as shown in [8]. For this reason, the effect of this 
approximation is essentially to replace the continuum of a interval of size 2L 
on the real line with a discrete set of — 1 points, Xk, uniformly spaced on 
this interval. 

Clearly these relations are easily generalized to functions of two or more vari- 
ables. Since the focus of this paper is on two dimensional membranes, we will 
briefly discuss how the LSF are used to discretize a region of the plane; the 
extension to higher dimensional spaces is straightforward. A function of two 
variables can be approximated in terms of (A^, — 1) x (A^ — 1) functions, cor- 
responding to the direct product of the — 1 and A^ — 1 LSF in the x and y 
axis: each term in this set corresponds to a speciflc point on a rectangular grid 
with spacings and hy (in this paper we use a square grid with = Ny = N 
and = Ly = L). 

Since {k, k') identifles a unique point on the grid, one can select this point 
using a single index 
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K = k' +- + {N -l)ik + --l 



(6) 



which can take the values 1 < K < {N — 1)^. This relation can be inverted to 



give 



K 



k' = K - N/2 - {N -1 



(7) 



K 



N-l + e. 



where [a] is the integer part of a real number a and e 0. 

To illustrate the collocation procedure we can consider the Schrodinger equa- 
tion in two dimensions: 



Hipnix, y) = [-A + V{x, y)] ipn{x, y) = E^ipnix, y) 



(9) 



using the convention of assuming a particle of mass m = 1/2 and setting 
h= 1. The Helmholtz equation, which describes the vibration of a membrane, 
is a special case of ([9]), corresponding to having V"(x, = inside the region 
B where the membrane lies and V{x^ ?/) = oo on the border OB and outside 
the membrane. 

The discretization of eq. do]) proceeds in a simple way using the properties 
discussed in eqs. Q and (|5j: 



Hkk' 



kk' 



Ckj^k'f + hjCk'jA + SkjSk'j'V{xk,yk' 



J2) 



(10) 



where {k,j, k',j') = —N/2 + 1, . . . , N/2 — 1. Notice that the potential part of 
the Hamiltonian is obtained by simply "collocating" the potential V{x,y) on 
the grid, an operation with a limited computational price. The result shown 
in (10) corresponds to the matrix element of the Hamiltonian operator H 



between two grid points, {k,k') and which can be selected using two 

integer values K and J, as shown in (|6|). 



Following this procedure the solution of the Schrodinger (Helmholtz) equation 
on the uniform grid generated by the LSF corresponds to the diagonalization 



of a {N — 1)^ X (A^ — 1)^ square matrix, whose elements are given by eq. (10). 
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2 Applications 



In this Section we apply our metliod to study tlie vibration of different non- 
homogeneous membranes. Our examples are a rectangular membrane with 
a linear and oscillatory density, a rectangular membrane with a piecewise 
constant density, a circular membrane with density p{x, y) = 1 + a/x^ + 
and a square membrane with a variable density which fluctuates randomly 
around a constant value. All the numerical calculations have been performed 
using Mathematica 6. 



2.1 A rectangular membrane with linearly varying density 



Our first example is taken from [T] and later studied by different authors 
[2]|3E] ; these authors have considered the Helmholtz equation over a rectangle 
of sides a and b, and with a density 



p{x,y) = l + a(^^ + ^^ . (11) 



Notice that the factor 1/2 appearing in the expression above derives from our 
convention of centering the rectangle in the origin, whereas the authors of [Tp] 
consider the regions x G (0, a) and y G (0, b). 

The Helmholtz equation for an inhomogeneous membrane is 



- Au{x, y) = uj'^p{x, y)u{x, y) (12) 

where u{x,y) is the transverse displacement and uj is the frequency of vibra- 
tion. 

As explained in [lOj, the collocation of the inhomogeneous Helmholtz equa- 
tion is straightforward, and in fact it does not require the calculation of any 



integral. The basic step is to rewrite eq. (12) into the equivalent form 



Au{x,y) = u^u{x,y) . (13) 



p{x,y) 



The operator O = —jj^Ty^ A is easily collocated on the uniform grid gener- 
ated by the LSF, and a matrix representation is obtained. To see how this is 
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10.24717742 


10 


12.64211438 


12.64210427 


10.62452059 


10.62409026 



Table 1 

Results for the first 10 frequencies for b/a = 1 and a = 0.1 (second and third 
columns) and to a = 1 (fourth and fifth columns). 

achieved we can limit ourselves to a one dimensional operator and make it act 
over a single LSF: 



\ cP 1 
-T^-r^Sk{h, iV, x) = - ^ — — c2^Sj(/i, iV, x)si{h, N, x) 
p[x) dx^ ^ p{xj) 

^-Et7^4'>.(/^,A^,^) • (14) 



The matrix representation of this operator over the grid can now be read 
explicitly from the expression above. It is important to notice that in general 
the matrix representation of O will not be symmetric, unless the membrane 
is homogeneous. From a computational point of view the diagonalization of 
symmetric matrices is typically faster than for non-symmetric matrices of 
equal dimension. 

In Table[l]we display the first 10 frequencies of the square membrane {h/a = 1) 
for a = 0.1 (second and third columns) and for a = 1 (fourth and fifth 
columns). The number of collocation points is determined by the parameter 
which is fixed to 10 (second and fourth columns) and to 12 (third and fifth 
columns). The comparison of the results corresponding to different gives 
us an information over the precision of the results: looking at the Table we 
see that typically the results agree at least in the first 5 digits, although we 
are working with a rather sparse grid. Also it should be remarked that the 
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method is providing a whole set of eigenvalues and eigenvectors, (A^ — 1)^ to 
be exact, whereas in other approaches each mode is studied separately. 

To allow a comparison with the results of ref. pl2]l3E] a calculation of the 
fundamental frequency of the rectangular membrane for different sizes of the 
membrane and different density profile is reported in Table [2j The numerical 
results have been obtained working with N = 12. These results can be com- 
pared with those of Table 1 and 2 of ref. [1] , of Table 1 of ref. p] , of Table 1 
of ref. [3] and of Table 10 of ref. |1] (in the last two references only the case 
a = 0.1 is studied). Comparing our results with those of Masad we have been 
able to confirm the observation of Laura et al., that the frequencies calcu- 
lated by Masad for a = 1 and b/a = 0.6,0.4,0.2 are incorrect. Actually, the 
results reported by Masad in these cases are just the second frequency (for 
b/a = 0.6,0.4) and the fourth frequency (for b/a = 0.2) of the corresponding 
rectangular membrane, which the author failed to identify as such. This ob- 
servation illustrates the advantage of working with a method which provides 
a tower of frequencies at the same time. 

Another great advantage of our method is the great rate of convergence which 
is typically observed as the number of grid points is increased. In the cases 
studied in ref. [TU] where the boundary conditions are enforced exactly (a circle 
and a circular waveguide) we observed that the leading non-constant behaviour 
of the eigenvalues for ^ 1 was 1/N'^. This analysis in now repeated here 
for the case of a = 1 and b/a = 0.2: the points in Figjl] correspond to the 
square of the fundamental frequency calculated using different grid sizes. The 
curves which decay with N correspond to fitting the numerical points with 
fr{N) = cq + Ci/N"^ with r = 3, 4, 5 respectively. The horizontal line is the limit 
value of fiiN) for oo, which corresponds to a; = 12.548431091 (notice 

that the result obtained for = 12 agrees in its first five digit with this result, 
whereas the analogous result of ref. [2] agrees only in three digits). The reader 
will certainly notice that f^iN) fits excellently the sets, thus confirming the 
observations made in ref. [TU]. A further important observation concerns the 
monotonic behaviour of the points in the figure: as observed already in ref. [TD] 
this method typically provides monotonic sequences of approximations which 
approach the exact value from above. 

In Fig. [2] we have plotted the first 200 frequencies of a square membrane with 
a = 0, 0.1, 1 (going from top to bottom), using a grid with A^ = 26. 



2.2 A rectangular membrane with oscillating density 

As a second example, we consider a rectangular membrane with density p{x) = 
1 + 0.1sin7r(a; -|- 1/2). This problem has been studied in ref. [SP] . In Table [3] 



7 



0/ a 


a = U.i 


= 1 


i 


A O O C! O O /I OO 'V 

4:.66o6o4:227 


3.610490268 


U.o 


A Qn71 fifiQ/l H 
4.yU ( ioDo4o 


4.UoiOiOoo 


0.6 


5.957896353 


4.942230449 


0.4 


8.252203964 


6.797302677 


0.2 


15.61334941 


12.54867967 



Table 2 

Results for the fundamental frequency of a rectangular membrane with density 
p{x) = 1 + a{x + 1/2) for a = 0.1 and a = 1, using N = 12. 
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Fig. 1. Square of the fundamental frequency of the rectangular membrane for 
b/a = 0.2 and a = 1 as a function of N. The dashed, solid and dotted lines cor- 
respond to fitting the numerical points with fr{N) = cq + ci/N'^ with r = 3,4,5 
respectively. The horizontal line is the limit value of fi^N) for N oo, correspond- 
ing to Ei=Lof = 157.4631229. 



we compare our results (LSF) with the results of ref. [3], for different sizes of 
the membrane. Our results, obtained with a grid corresponding to = 20, 
agree with those of ref. [3j . In Table |4] we compare the results for the first 
10 frequency of the square membrane with this density with those reported 
in ref. Although our first 7 results agree with those of ref. we have 
noticed that the remaining results do not agree. By looking at the table the 
reader will notice that the disagreement is caused by the fact that ref. |1] has 
missed a frequency (such an error cannot take place with our method, since the 
diagonalization of the Hamiltonian matrix automatically provides the lowest 
part of the spectrum). 
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(I) 




Fig. 2. First 200 frequencies of a square membrane with density p{x) = l+a(x+l/2) 
for a = 0, 0.1, 1 (from top to bottom) obtained using N = 26. 
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LSF 


Ref. [3] 
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Table 3 

Results for the fundamental frequency of the rectangular membrane with density 
p{x) = l + asin7r(x + l/2) for a = 0.1, using the LSF with = 20 (second column). 
The third column are the results of Ref. [3] . 
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Fig. 3. Fundamental frequency of the rectangular membrane with discontinous den- 
sity. The solid lines are the fits f{N) = cq+ci/N over the three different monotonous 
sequences of numerical approximations. 
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Table 4 

First 11 frequencies of a square membrane with density p{x) = 1 + 0.1 sin tt{x + 1 /2) 
using the LSF with = 20 (second column). The third column are the results of 
Ref.t4j. 

2.3 A rectangular membrane with discontinous density 

Our next example is taken from ref . [5|6] : it is a rectangular membrane of sizes 
a = 1 and h = 1.8. The membrane is divided into two regions by the line 



y + b/2 = 0.3{x + a/2) + 0.7 . (15) 

The upper region has a density which is twice as big as the density of the lower 
region. The collocation procedure is the same as in the previous example. In 
Fig. [3] we have plotted the fundamental frequency of the membrane calculated 
with different grid sizes, i.e. different A^. In this case we clearly observe an 
oscillation of the numerical value: as discussed in ref. [TU] this behaviour is 
typical when the grid does not cross the boundary. Although in the present 
case the Dirichlet boundary conditions are enforced exactly on the border 
of the rectangle, the region of discontinuity is not sampled optimally by the 
grid, which causes the oscillation. Nonetheless, we may observe that the set 
of numerical values can be divided into three distinct and equally spaced sets, 
each of which can be fitted quite well with a behaviour f{N) = cq + ci/N (the 



We use the convention of centering the membrane in the origin of the coordinate 
axes. 
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Table 5 

First 10 frequencies of a rectangular membrane with discontinous density for differ- 
ent grid sizes. 

solid curves in the plot). Going from top to bottom, the curves correspond to 
the fits: 

/(iV) = 2.76810 + 0.0975751/A^ (16) 
/(AT) = 2.76778 + 0.0395891/A^ (17) 
f {N) = 2.76785 - 0.0304169/A^ . (18) 



2.4 A circular membrane with density p{x, y) = 1 + + 

Another interesting example is taken from ref. [4j, where an inhomogeneous 
circular membrane with density p(x, ?/) = ! + ^/x^~+^ is considered. Ta- 
ble 8 of ref.[l] contains the first 5 eigenvalues. We have applied our method 
to this problem, using grids of different sizes, with going from 10 to 30. 
Studying the A^ dependence of the eigenvalues, we have seen that these de- 
crease monotonically and that the leading non-constant dependence on A^ for 
A^ — >■ cxo is A^~^ (see Figjs]). In Table |6] we report the first 10 frequencies of 
a circular membrane with density p{x, y) = 1 + a/x^ + for different grid 
sizes (A^ = 26,28,30). Notice that the results already agree in their first 4 
digits. A more precise result is then obtained by performing an extrapolation 
of the numerical results for grids going from A^ = 12 to A^ = 30. Notice the 
good agreement with the results of (although we believe that our results 
are more precise), and that, as expected, some frequencies are degenerate (the 
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Fig. 4. Density of the square inhomogeneous membrane obtained by conformal map- 
ping of the circular homogeneous membrane with density p{x, y) = 1 + \/ + y^- 
The sohd hues mark the level 1,2 and 3. 
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Fig. 5. Square of the fundamental frequency of the circular membrane with density 
p{x^ y) = 1 + as a function of N . The dashed, solid and dotted lines 

correspond to fitting the numerical points with fr{N) = cq + ci/N"' with r = 3, 4, 5 
respectively. 

degeneracy of the frequencies is not discussed in [1]). 
2.5 A square membrane with a random density 

Our last example is a square membrane with a variable density which fluc- 
tuates randomly around the value Pq = 1. We have generated random values 
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Table 6 

First 10 frequencies of a circular membrane with density p{x, y) = 1 + \/ + y"^ 
for different grid sizes. The fourth column reports the frequencies obtained by ex- 
trapolating the numerical results for grids going from = 12 to A'^ = 30. The last 
column reports the results of ref. [3] . 

for the density on a uniform square grid with 81 points (corresponding to 
A^o = 10 in our notation); at these points the value of the density has been 
chosen according to the formula: 

p{xk,yj) = po + Sp Qkj (19) 

where po = 1 and 6p = 1/2. q^j is a random number distributed uniformly 
between —1/2 and 1/2. The density over all the square has then been obtained 
interpolating with the LSF: 

p{x,y) = Cj2pi^k,yj) Sk{h,No,x)sj{h,No,y) , (20) 

where A^o = 10, as previously mentioned. C is a normalization constant which 
constrains the total mass of the membrane to be equal to the mass carried by 
the homogeneous membrane with density pQ. 

In the left panel of Fig|6]we plot the density of the membrane, while in the right 
plot we plot the fundamental mode. We have performed our calculation using a 
grid with = 40 (i.e. with a total of 1521 modes): upon diagonalization of the 
matrix obtained using the collocation procedure we have obtained numerical 
estimates for the lowest part of the spectrum of the random membrane. These 
results may be compared with those of the homogeneous square membrane 



13 



and with the asymptotic behavior predicted by Weyl's law [TT] : 



-'-^^jf-^, (21) 

A being the area and L the perimeter of the membrane. In Fig. [7] we plot the 
quantity A„ = En — E^^^^ for the homogeneous square membrane (squares) 
and the square membrane with randomly oscillating density (circles). We have 
limited the plot to the first 200 modes. Notice that in both cases A oscillates 
around 0, although the oscillation are smaller for the random membrane. 

We have also fitted the first 400 modes with the functional form given by 
Weyl's law obtaining 



En ^ 3. 14546n + 3.45835 x/n (22) 
for the homogeneous membrane and 



En ^ 3.14445n + 3.44674 Vn (23) 
for the random membrane. These values should be compared with the one 



given by eq. (21 ): 



^ 3.14159n + 3.54491v^ . (24) 



3 Conclusions 



In this paper we have introduced a new method to solve the Helmholtz equa- 
tion for non-homogeneous membranes. This method uses the Little Sine Func- 
tions introduced in refs. pi9] to obtain a representation of the Helmholtz 
equation on an uniform grid. The problem thus reduces to diagonalizing a 
(A^ — 1)^ X (A^ — 1)^ square matrix, — 1 being the number of collocation 
points in each direction. We have tested the method on several examples taken 
from the literature. The application of the method is straightforward and it 
provides quite accurate results even for grids with moderate values of A^. The 
readers interested to looking to more examples of application of this method 
should check the site 

\protect\vruleuWidthOpt\protect\href {http : //f e j er . ucol . mx/paolo/ drum}-[http : //f e j er . uc 
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Fig. 6. Left PanehDensity of the square membrane with random density (Set 1). 
Right Panel: Fundamental mode of the square membrane with random density (us- 
ing Set 1). We have used the LSF on a grid with = 40. 




Fig. 7. An = En — En '^^ for the homogeneous square membrane (squares) and the 
square membrane with randomly oscillating density (circles). 

where images of modes of vibration of membranes of different shapes can be 
found. 



References 

[1] J.A.Masad, Free vibrations of non-homogeneous rectangular membrane, Journal 
of Sound and Vibration 195, 674-678 (1996) 

[2] P.A.A.Laura, R.E.Rossi and R.H.Gutierrez, The fundamental frequency of non- 



15 



homogeneous rectangular membranes, Journal of Sound and Vibration 204, 373- 
376 (1997) 

[3] S.H.Ho and C.K.Chen, Free vibration analysis of non-homogeneous rectangular 
membranes using a hybrid method. Journal of Sound and Vibration 233, 547-555 
(2000) 

[4] S.Yu. Reutskiy, The methods of external and internal excitation for problems 
of free vibrations of non- homogeneous membranes. Engineering Analysis with 
Boundary Elements 31, 906-918 (2007) 

[5] S.W.Kang and J.M.Lee, Free vibration analysis of composite rectangular 
membranes with an oblique interface. Journal of Sound and Vibration 251, 505- 
517 (2002) 

[6] C.P.Filipich and M.B.Rosales, Vibration of non-homogeneous rectangular 
membranes with arbitrary interfaces. Journal of Sound and Vibration 305, 582- 
595 (2007) 

[7] P. Amore, A variational sine collocation method for strong-coupling problems. 
Journal of Physics A 39, L349-L355 (2006) 

[8] P.Amore, M.Cervantes and F.M.Fernandez, Variational collocation on finite 
intervals, J.Phys.A 40, 13047-13062 (2007) 

[9] P.Amore, Alternative representation of nonlocal operators and path integrals, 
Phys.Rev.A 75, 032111 (2007) 

[10] P. Amore, Solving the Helmholtz equation for membranes of arbitrary shape, 
sent to Journal of Physics A (2008) 

[11] J.R.Kuttler and V.G.Sigillito, Eigenvalues of the Laplacian in two dimensions, 
Siam Review 26, 163-193 (1984) 



16 



